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Abstract 

This paper presents a simple ball-and-spring model for the propagation of small ampli- 
tude vibrations in a granular material. In this model, the positional disorder in the sample 
is ignored and the particles are placed on the vertices of a square lattice. The inter-particle 
forces are modeled as linear springs, with the only disorder in the system coming from a 
random distribution of spring constants. Despite its apparent simplicity, this model is able 
to reproduce the complex frequency response seen in measurements of sound propagation 
in a granular system. In order to understand this behavior, the role of the resonance modes 
of the system is investigated. Finally, this simple model is generalized to include relaxation 
behavior in the force network - a behavior which is also seen in real granular materials. 
This model gives quantitative agreement with experimental observations of relaxation. 



PACS numbers: 63.50,+x, 05.40.+j, 46.10,+z, 43.40,+s 



Jan. 25, 1993 



1. Introduction 

Understanding the properties of granular materials present a difficult challenge to 
condensed matter science 0. A fundamental property of any material is its response to 
small amplitude vibrations. For the case of solids, the behavior of the system is character- 
ized by the phonon spectrum. These elementary excitations for perfect crystals are well 
understood, but for disordered systems, the situation is not so clear. When disorder is 
present, localization of these excitation can occur ||, and while much progress has been 
made in understanding the phenomena of localization, there are still many fundamental 
properties that remain a mystery ||. 

While there have been many studies of the localized modes in disordered systems 
J§-|(|, there has been less attention paid to the response of these systems to a periodic 
driving force 0. Generally speaking, for a uniform periodically-driven linear system in 
the long time limit, the frequency of oscillation is that of the forcing ||. As the frequency 
of the driving goes through a resonant frequency, the amplitude of the motion is increased, 
and the oscillations takes on the character of the corresponding normal mode. Does this 
picture remain for the case of a disordered granular material which has complicated non- 
linear inter-particle interactions? In these systems, there is known to be a network of 
contact forces with a complex structure. The phenomena of arching suggest that this 
system is supported by a small fraction of grains which have a high concentration of stress, 
while most of the particles remain in loose contact 0. The response of this system to 
driving may then depend wholly on the structure of this network where the stress is highly 
concentrated. 

In this paper, I discuss a model system of balls and springs that can reproduce many 
of the vibrational properties of a real granular system. The numerical results that come 
from this model will be compared with the results from a beautiful set of experiments 
performed by Liu and Nagel in which they measured the response of a single "grain of 



sand" to an external driving force |I(J . I describe briefly these experiments and summarize 
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the important results. The experimental apparatus consisted of a rigid box filled with glass 
bead of diameter d = .5 cm. The box itself was 28 cm x 28 cm and was filled with beads 
up to heights of 8 to 15 cm. The upper layer of the bead pack was a free surface. At one 
end of the box was placed an aluminum disk that was attached to an external speaker by a 
rigid rod. The speaker could be driven with varying frequencies and amplitudes. Inside the 
bead pack was an accelerometer which was roughly the size of one bead. This device was 
sensitive only to horizontal accelerations, and was unaffected by sound waves propagating 
through the surrounding air. Since this detector was comparable to the size of a bead, 
it is effectively measuring the motion of a single particle under the action of the driving 
vibration. The motion of the aluminum disk was also monitored with an accelerometer 
attached directly to it. 

The disk was driven with an acceleration of the form A s sin(wt) and the detector was 
found to oscillate as Adit) sin(o;i + </>(£)), where Ad(t) and cf)(t) were the detector amplitude 
and phase shift, respectively. Both of these functions were found to vary slowly with time. 
Figure |l|(a) shows Ad(t) with a driving frequency uj = 4 kHz, and A s = lAg where g is the 
acceleration of gravity. Note that the time scale for the changes in Adit) is much longer 
than the time scale for the oscillation itself. The power spectrum (shown in fig. [j](b)) of 
this time series shows a power law region with exponent ~ 2 at frequencies from 10 -5 Hz 
to 10" 1 Hz. 

On time scales of the order of a few oscillations, however, the amplitude at a given 
driving frequency is fixed. This inspired a second set of measurements using the same 
experimental set-up. With a fixed acceleration for the driving force, the amplitude of the 
detector is measured for various value of uj. Let the response function, rj(u) = Ad/A s for a 
driving frequency, uj. This ratio is plotted as a function of uj in fig. ^. The two curves shown 
were done in consecutive sweeps through the frequency range with the measurements being 
separated by a short interval in time. The second curve is displaced vertically downward 
so that the two curves can be distinguished. If the sample was disturbed between the 



4 

measurements, the response curve looked very different for the second scan of the frequency 
range. Thus, while the data looks "noisy," the response is actually characteristic of the 
particular contact network of the sample. Thus, fig. || gives information about a static force 
network in the bead pack. Figure [I[ on the other hand, results from the time evolution of 
this network. 

In the rest of this work, I present the details of a ball and spring model for the bead 
pack and present numerical results for a direct comparison between the model system and 
the experimental results. I also study the role that the system's normal modes play in the 
frequency response. The rest of the paper will be organized as follows. Section 2 contains 
a discussion and motivation for a ball and spring model for this system. In section 3, I 
will consider the continuum limit of this model and examine the solution to the differential 
equations which describe a homogeneous system. Section 4 presents the numerical solution 
of the discreet model, and compares these results with the frequency response seen in the 
experiment. In section 5, I present a toy model for the slow amplitude modulations which 
are a result of the relaxation of the force network. Finally, in section 6, I summarize 
and discuss extensions to this simple model which must be made to capture more of the 
physical characteristics of the real system. 

2. Ball and Spring Model 

Consider the bead pack at equilibrium. In this experiment, the only important forces 
acting on each bead are the contact forces between adjacent particles and gravity. The 
collection of contact forces among neighboring beads will be referred to as the contact force 
network, and I sometimes refer to the force between two connected beads as a "bond." 

The fundamental unit of this network is the bond between two spheres. If two spherical 
beads are compressed together with a contact force F c , then under the action of this force 
the distance between the centers of the two spheres is decreased by an amount h. For the 
case of linear elasticity and perfectly smooth spheres, these two quantities are related by 

F c = kh 3 / 2 , (2.1) 



where k is a proportionality constant depending on the radii of the two spheres and their 



material properties |TT| . For the case of an almost mono-disperse collection of beads, this 
constant k can be treated as the same for all contacts. 

Consider now an additional force, SF, acting on these two balls. The spheres will 



move together by an additional amount 5h. Under this new force, equation ( [2.1| ) becomes 

F C + SF = k(h + Sh) 3/2 . (2.2) 

For the case of small Sh, I can expand the expression in parenthesis and making use of Eq. 
(0), I find 

SF = ^—5h = K h Sh, (2.3) 

where Kh is a constant which depends on the equilibrium stress on each bond in the 
network. Equation (|2.3| ) shows that for a sufficiently small displacement, 5h, the bond 
between a pair of connected beads acts like a simple spring. 

Unfortunately, ( |2.3|) is not true for all bonds. There may exist some beads which in 
equilibrium are very close to each other, but are not in contact. However, under the action 
of 5F these beads do come into contact. The experimental results indicate that this is not 
an important effect in understanding the response of the system. These contacts produce a 
force that only acts during part of the oscillation cycle. If this effect were important, then 
the response of the accelerometer would not be sinusoidal, but would take on some more 
complicated wave form. Thus, in order to understand the experimental results described 
above, I will not consider these "sometimes" bonds. 



Finally, one might object to the use of a Eq. (|2.1|) on the basis that it ignores the 
effect of the surface roughness of the beads. However, the essential idea that for small 
oscillations the contact will behave like a spring does not depend on the fact that the 
contact be of the Hertzian type. All that is required is an analytic form for the behavior 
near equilibrium in order to do an expansion for small 8h. The dependence then of Kh 
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on the equilibrium configuration may then become more complicated; however, I will not 
need the analytic form for in this study. 

In terms of the response to small oscillations, the system can then be viewed as a 
collection of balls connected by springs, with the spring constants determined in some way 
by the force network of the equilibrium system. When the amplitude of a ball changes as a 
function of time, this indicates that the equilibrium configuration of the system has changed 
- the network has relaxed. The experimental results then raise some very interesting 
questions: What information about the force network is reflected in figs. [I] and [2]? How 
does the network relax in time? What is the role of the geometrical disorder in the system? 
These are the questions that I address in this study. 

I focus on an an extremely simple model for the system. I consider a network embedded 
in two spatial dimensions. I ignore the positional disorder, and consider the network 
topology to be that of a square lattice. I choose a velocity dependent dissipation with a 
uniform damping constant. The only disorder present in the system is in the distribution 
of spring constants, K^. This paper will compare the response of this model system to the 
response observed in the experiments. 

I now present the mathematical formulation of this ball and spring model. Consider 
the beads, confined to a two dimensional square lattice as in fig. |3|. I let the driving force 
be a sinusoidal motion of one entire wall. Because energy is constantly being pumped into 
this system by the vibration, there must be some method for dissipation in the system. I 
assume a damping force that is linearly dependent on the velocity difference between two 
adjacent beads. The final simplification comes from an assumption of the decoupling of the 
motion in the x and y directions. This will be true in the limit of very small oscillations. 

In this model, the displacement u(p) of a bead at position p = (x, y) is given by the 
differential equation 
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where m is the mass of a ball, p' is a nearest neighbors to point p, (3 is the damping 
constant, and Kh(p,p') = Kh(p',p) is the spring constant between p and p' . 

For boundary conditions, I take that that u(p) = on all walls that do not supply a 
driving force. For the oscillating wall, I have the form u = exp(iut), where, without loss 
of generality, I have taken the amplitude of the wall's vibration to be unity. 

It is possible to take Eq. ( |2.4j ) as the starting point for numerical investigation by 
merely integrating this equation to find out the response of each ball as a function of 
frequency. However, this does not take advantage of the fact that for small amplitude 
vibrations, the harmonic driving force produces a harmonic motion of the particles in the 
bead pack. Thus, I know that the solution of interest takes the form u(p) = A(p) exp(zwt), 
where A(p) is the complex amplitude for the vibration at point p. Combining this with 
Eq. 0, I find 

mu 2 A{p) + J2(K h (p,p') + (A(p') - A(p)) = 0. (2.5) 
p> 

The boundary conditions for this equation is that A(p) = on the three fixed boundaries, 
and that A(p) = 1.0 on the oscillating wall. This equation described the steady state 
response of this ball and spring model. To find the frequency response for a given system, 
it is necessary to know the values of Kh(p,p'). 

3. The Continuum System 

If I consider the case of infinitesimally small particles (i.e. the continuum limit), it is 
possible to rewrite Eq. (|2.5|) in terms of a differential equation. The resulting equation can 
be solved exactly for the case of a uniform spring constant. Letting Kh(p,p') = K for 
all p and p' , I write A(p) as A(x, y) and consider a box with x and y dimensions L x and 
width L y , respectively. For such a system, ( |2.5|) becomes 

muj 2 A(x, y) + (K + iu;(3)V 2 A(x : y) = 0. (3.1) 



The boundary conditions for this equation are A(L x ,y) = 1.0, and A(x,0) = A(x,L y ) = 
A(0, y) = 0.0. The solution to this equation, with these boundary conditions can be 
obtained by standard separation of variable techniques, and I find 

K J *^ (2n + 1) V sinhK) y ' 1 J 

where a n is defined by 



a ^ = (2 „ + 1 )V/ i; -^. (3.3) 

It is very simple to sum this series numerically, and find the amplitude as a function 
of frequency for any point in the box. All that is necessary is to choose values for the 
parameters m, Kq, /3, L x and L y . In the solutions shown here, I choose m = 1, Kq = 1, 
L x = L y = 1. Figure § shows the phase and frequency response at (x,y) = (.47, .52) 
for two different values of f3. In fig. f|(a), I show the results for (3 = 0. There is a 
set of frequencies where the amplitude at (x, y) is peaked, and is much greater than the 
amplitude of the driving wall. Additionally, there are frequencies at which the amplitude 
jumps discontinuously. The frequencies, u p , at which these peaks and discontinuities occur 
are well described by the expression 

u 2 p = kn 2 (n 2 l + (2n 2 + l) 2 ), (3.4) 

where n\ > and n<i > are both integers. The frequencies defined by (|3.4| ) are shown in 
fig. ||(a) by the pluses located on the x-axis. 

The importance of these frequencies can be understood in two ways. The first comes 
from a direct examination of ( |3.2j ) and (373). At uj = uj p , the value of a n2 = riini. For 



this value of a n , the function sinh(a n ) vanishes, and the contribution for that the n = n 2 
term in the series ( |3.2| ) goes to infinity. 
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However, the frequencies of (|3.4|) also correspond to natural frequencies of the unforced 
system. Let (3 = and A(L X , 1) = 0. In this case, Eq. ( j3.1| ) becomes an eigenvalue problem, 
and the eigenfunctions, fi(x,y), are 

fi(x, y) = sin(2q 1 irx/L x ) sin(2q 2 7iy / 'L y ), (3.5) 

with the corresponding eigenvalues, uf, given by 

^ = kn 2 (ql + q 2 2 ), (3.6) 

where again qi > and q 2 > are both integers. Not all eigenfrequencies given by Eq. 
( |3.6|) are manifested in the system. This is due to the symmetry of the driving force - 
it is symmetric about a horizontal axis at the center of the box. This excludes all of the 
eigenmodes which are asymmetric about this axis (i.e. have an even value for q 2 ). 

Figure |4|(b) shows the response for a finite value of the damping coefficient, (3 = .003. 
The frequencies defined by ( |3.4| ) are again shown as pluses along the x-axis. There are 
two distinct differences between this response curve and fig. J§(a). The first is that there 
is an overall decrease in response as the frequency increases. This reflects the fact that 
the damping force is proportional to u, and thus the damping force gets larger at higher 
frequencies. By looking at a larger range in frequency, it is clear that this trend is nothing 
as straight-forward simple exponential decay. 

The second interesting feature is the set of peaks and valleys which are superimposed 
on this decreasing response curve. The lowest frequency excitations are still identifiable 
as distinct entities with peaks clearly occurring at the resonance frequencies. The higher 
frequency modes, on the other hand, can not be distinguished in the response curve. 
The finite value of the damping term causes the resonance peaks to be greatly decreased 
in magnitude, and to be much broader in terms of their frequency response. At higher 
frequencies, many broad peaks are superposed to get the measured response of the system. 
At some frequencies the result is a slightly higher than average response, and at other 
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frequencies, it is slightly lower. This causes the apparent small peaks and valleys at high 
frequencies. At even higher frequencies, the effects of the damping become dominant, and 
the curve becomes smooth. 

This is indeed reminiscent of the frequency response in the experimental system. 
There is a general trend of a decreasing amplitude with increasing frequency, with peaks 
and valleys superimposed on this decreasing curve. In the experiments, these oscillations 
are much more striking than in this uniform system, but the behavior is the same. 



4. The discreet model 



I now study Eq. ( |2.5| ) numerically, looking at the properties of the solution for a par- 
ticular distribution for the values of Kh(p,p'). The method that I use to find the solutions 
of equation ( |2.5| ) is the Bi-Conjugate Gradient Method fl2]] . This method determines 
solutions to a matrix equation of the form 

Mx = 6, (4.1) 

by finding the minimum of the function g(x) — ■ M ■ x + b ■ x. In order to transform 
Eq. (|2.5| ) into the form of (O), it is necessary to write (|2.5| ) as an equation for its real 



and imaginary parts. As a result, the matrix M is not symmetric, and this is why the bi- 
conjugate, rather than the conjugate, gradient method must be used. In my simulations, 
the amplitude for each point is accurate to 1 part in 10000. 

It is, of course, possible to simply invert the matrix A to find the solution. However, 
this does not take advantage of the sparseness of the matrix. All of the calculations for the 
results shown in this paper were done using a Sun SPARCstation IPX using approximately 
150 hours of CPU time. 

All that remains is too choose a distribution for the set of {Kh(p,p')}. For simplicity, 
I take a random distribution of values on the interval [0, K max ] . The results that I show 
are for a square array of balls of size 20 x 20. Figure [5| shows the spring configuration 
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for the array that has been used for the results presented below. Black indicates that 
Kh = K max , while white means Kh = 0, with a linear grayscale for intermediate values. 

I choose m = 1, K max = .25, (3 = .003 and measure the response in the frequency 
range from 10 = (.01, 1.3). I shall show below how to assign physical units to these values. 
Figure ^| shows the amplitude for a ball near the center of the pack of the springs (the 
precise position of the ball is indicated in fig. |5|). 

There are three distinct regimes of behavior. At very low frequencies (fig. §|(a)), the 
behavior is very similar to the response seen at low frequencies in the continuum case. 
There are frequencies at which peaks occur, and these peaks show little overlap. 

In the intermediate frequency range (fig. 0(b)), the response of the system is quali- 
tatively the same as that shown in the experimental results of fig. |^. There is an overall 
decrease of the response, and an irregular set of peaks and valleys. 

For u > 1 (fig. H(c)), the response becomes very smooth and decreases very quickly. 
The change in response at u ~ 1 can be understood intuitively. Consider a single ball 
attached to four springs all with Kh = K max . Such a ball will have its fundamental 
excitation at a frequency 



This then is the "natural" frequency (known as the Einstein Frequency ||13|| ) for a ball 
which, by chance, happens to be surrounded by springs all of which have Kh ~ K max . 
Equation (|4.2|) then defines the highest frequency at which any bead will naturally oscillate. 
Any higher frequency oscillations would require every neighbor to vibrate in opposition 
with many of its neighboring particles. Below it is shown that this results in a distinctly 
different character to the system's oscillations. 

For the values of K max and m given above, iv max = 1.0. For other simulations, 
with different values for K maxi I found that this change in the frequency response of the 
system occurred always at the frequency given by Eq. (|4.2| ). This frequency determines 
the frequency scale for the system, and I will consider my frequency as measured in units 




(4.2) 
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of Umax. Experimentally, u) max might be measured by finding the frequency at which this 
rapid decrease in the response occurs. It is also possible to estimate a lower limit for u> max 
in the experimental system by assuming the Hertzian contact law and that contact force 
between any two beads is not significantly greater than the weight of a single bead. Making 



use of some of the information about the beads given in ref. [10], and Eqs. (|2.1| ), ( |2.3| ) and 
(0), I find that 

oJmax > 45kHz, (4.3) 

which is above the observed range in the experiment. 

Consider now the phase of a given ball in the array. Figure [F] is a plot of the phase 
versus frequency for the same ball as was used in fig. |6|. Again there is distinctly different 
behavior for u > oj max . There is also a linear regime in the phase for .3u max < u < 
Umax- A linear relationship between the observed phase and the driving frequency is seen 
experimentally [TTJ]]. A simple calculation shows that for a constant phase velocity, v^, the 



relationship between the frequency of the oscillation, and the observed phase at a point 
should be linear. If the slope of the phase response is fi, then the phase velocity is given 
by 

= X/n, (4.4) 

where X is the distance from the source of the driving to the point of observation. In the 
data shown here, the value of X is 12 particle diameters. In the experimental system, the 
particle diameter is .5cm. Using this value and the slope from fig. [7], I find that 

v<j, « .03u max , (4.5) 

where uj max is measured in sec -1 , and the velocity is given in units of cm/sec. In the 
experimental system, = 57m/sec. Matching to that velocity implies that 



UJ 7 



200,000sec" 1 , (4.6) 
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which is clearly beyond the experimental regime. With this value of uj max , fig. |6| indicates 
that the typical Au over which the amplitude changes appreciably is ~ 20 kHz. This 
compares with the experimental value which is quoted as Auj « 3 kHz. 

The parameters m and K max determine the value of iv max through (|4.2|) . The only 
parameter whose effect I have not discussed is the damping constant (3. The value observed 
for \i depends sensitively on the value chosen for (3. If a smaller value for (3 is chosen, the 
value of \i is increased. Additionally, a smaller f3 results in many more small peaks in 
the frequency response data, but the sharp decreases are much shallower. If, on the other 
hand, the value of f3 is increased, the sharp decreases in the response are much deeper, 
but they occur with a much larger spacing in frequency. The effects on the peaks can be 
understood my considering the results from the normal mode analysis below. 

Because the normal modes proved to be so important in the response of the continuum 
system, I consider their manifestation in the discreet case. I again use the configuration of 
fig. |5|. The normal modes are calculated for the undamped system with the zero amplitude 
boundary conditions at all walls. The numerical calculation of the eigenmodes is done 
using the appropriate subroutines from LINPAK. Figure |8| shows three eigenmodes for 
the system. The first one is the second lowest frequency mode of the system, and shows 
little trace of the local non-uniformities of the network. The second mode is in the middle of 
the spectrum at u> = .6053u; max . The picture shows some localization of the oscillation, but 
not complete. The final mode shown has ui = 1.000uJm ax , and there is striking localization 
of this mode. This is true of all the modes at high frequency, and this is the reason for the 
significant change in the response function at driving frequencies with ui > uj max . 

These normal mode frequencies are shown as pluses in figs. ^ thru [7[ Notice, as in the 
continuum case, that at low frequencies the sparseness of the normal mode allows them 
to be observed as distinct excitations of the system. At higher frequencies, however, the 
modes overlap and there can be no one-to-one identification between the structures in the 
response and a given normal mode. 
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Note also that the density of modes becomes very low for frequencies greater than 
u « uJm ax . From the results for low frequencies, this would seem to suggest that these 
modes should be seen distinctly in the response curve. However, because these modes 
are highly localized, they will only be seen in the response if the detector happens to be 
inside the region of space where the mode has a non-zero value. In addition, these highly 
localized modes must have some overlap with the boundary or they cannot be significantly 
excited by the driving force. 

The major difference in this non-uniform system is that not all modes will contribute 
equally to the system's response. A high frequency mode which has its amplitude localized 
far from the driving wall is much less important than a mode which is localized near 
the driving wall. Thus, not only are the functional forms of the normal modes more 
complicated, it is also important to know which modes make the dominant contributions 
at each frequency. 

In order to measure the contribution of each mode to the system response, I use the 
set of normal modes as a set of basis vectors to describe the amplitudes of the balls at 
each driving frequency. At a given driving frequency u, the amplitude of the oscillations 
is given by A(x, y). A(x, y) can be written as a superposition of the eigenmodes, 

NxN 

A(x,y) = c ih(x,y), (4.7) 

i=0 

where c$ is given by 

N N 

C^^A^fcj/^fc). (4.8) 

j=0 k=0 

Thus, by plotting the \ci\ as a function of frequency, I can see when each mode is con- 
tributing to the behavior of the driven system. 

With the modes calculated numerically, it is a simple matter to calculate the values of 
Ci for each driving frequency. Figure |^ shows the values of the \ci\ versus uji at a driving 
frequency ui = .6u)max- While it is clear that the most important contributions come from 
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the modes that have Ui ~ u, the structure of the curve is complicated. It is also clear that 
the level of a modes excitation is not simply determined by the value of \u — 

There are, however, some general patterns to these Cj's. Figure [10] shows the fre- 
quency of the mode which has the largest value of \ci\ at each driving frequency. This is 
a linear function with slope 1, and indicates that the most important mode is one which 
has its resonant frequency near the driving frequency. Again, it is clear that the behavior 
sharply changes at u ~ oj max . 

Finally, I consider the c^'s for several eigenmodes which have their frequencies very 
close together. Shown in fig. [11] are the values of the \ci\ as a function of frequency for three 
modes with a resonance frequency near uj = .6w max . Each curve shows a clear resonance 
behavior when the driving frequency approaches its natural frequency. However, it is clear 
that each mode is excited to a different degree by the driving force. 

As was stated above, it is possible to understand the changes in the response curve 
as the value (3 is changed by considering the normal mode. When the value of (3 is 
increased, each of these modes will have a much sharper peak in its when the driving 
frequency approaches the mode's natural frequency. Thus, any mode which has a non-zero 
contribution to the amplitude at the detector is more likely to make its presence felt. Thus, 
there are many more smaller peaks in the measured response. 



5. Toy Model for Relaxation 



Because ( |2.5|) is a steady state solution for fixed u, there will be only one amplitude 
measured for all time. I have allowed for no relaxation of the lattice which is the source of 
the behavior in fig. |]. As a toy model for this relaxation process, I consider the amplitude 
at one point, while randomly changing one bond per unit of time, r. This bond is assigned 
a new random value between and K max . Having changed one of the values of K^, I then 



find the new solution of (|2.5| ). This new solution has a different value for the observed 
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amplitude, and by plotting this as a function of elapsed time, I generate a time sequence. 
I then compare the power spectra of this series to the experimental power spectra. 

I begin the simulations with the same bonds configuration as shown in fig. |5|, and 
measure the amplitude of the same ball as I used for the frequency data. The system 
is driven with an oscillation of frequency 10 = .6c<; max , with a value of (3 = .01. The 
amplitude of the vibration as a function of time is shown in fig. |l2|(a). Using a much 



longer time series, I can calculate the power spectrum, S(f), for this time series jig] , with 
the results shown in figure fig. |l2](b). There is a clear power law behavior in the power 
spectrum. If the power spectrum is written as S(f) ~ / _ , then I find that tjj ~ 1.9 over a 
region of one and a half decades. This is fair agreement with the experiments, which find 
values for ijj between 2 and 2.2. 

It is not surprising that ijj ~ 2 in both the toy model and the real system for the 
following reason ||10|| . Imagine that at one time the amplitude of the detector has some 
value Bq and then one of the bonds is changed. This will either raise or lower the observed 
amplitude by ABq, and I might expect, in general, that the change might be equally likely 
to cause an increase as a decrease. Then another bond is changed, and the amplitude 
changes by some amount ± ABi . In this picture it is clear that the amplitude is executing 
a random walk. Such a time trace is known to have a power spectrum with an exponent 
ip = 2. Thus, the power law behavior is not surprising, and the exponent near 2 seems to 
be intuitive. 

There are, however, several assumptions in this argument which are unjustified. Are 
B and AB actually uncorrelated? Are successive values for AB uncorrelated? In order to 
understand the experimentally observed results, a more physically motivated process for 
the network relaxation is necessary. 



6. Conclusions and Discussion 

The principal objective of this paper is to show that the vibrational properties of a 
granular material can be reproduced with a simple linear ball and spring model. I have 
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been able to reproduce qualitatively the frequency response of a single grain to a sinusoidal 
driving frequency. There is as yet no way to quantitatively compare with experiments, 
because there is no well-defined method for characterizing the statistical properties of the 
"noisy" response curve. The normal mode analysis indicates that the peaks in the spectrum 
cannot be associated on a one to one basis with a resonance of the system; however, it is 
clear that the behavior at a given frequency is a results of a superposition of the normal 
modes which have their resonance frequencies in the regime of the driving frequency. The 
toy model for the relaxation of the lattice does reproduce a power law in the behavior in the 
power spectrum. The fact that the model is simple, and yet produces the experimentally 
observed behavior, suggests that it may be worth studying in its own right. 

However, there is at least one aspects of the real system that can not be studied 
with this simple ball and spring model - the experiments shows an extreme sensitivity to 
thermal fluctuations. In additional measurements, Liu and Nagel place within the sample 
a small heater the size of a single bead, and then run a heat pulse through it [[yj. The 



pulse changes the temperature of single bead by approximately 0.8 K, and they find that 
the response of the distant accelerometer changes almost instantaneously by 25%. The 
effect of the change in temperature is to change, by thermal expansion, the size of the 
beads in a region near the heater. The change in radius of the bead is estimated at 1000 
A. In order to understand these temperature effects, a much more physically motivated 
model is needed for the granular system. 

The major defect with the current model is that it does not allow for a general topology 
for the force network. Also, the spring constants are clearly determined by the different 
equilibrium forces on the beads. However, there is no requirement that these forces sum to 
zero, as they should in equilibrium. This problem can only be solved by choosing a specific 
form for the force law between touching spheres. Once these two problems are solved, it 
would be straight forward to study the effects of local thermal fluctuations. 
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Figure Captions 

Fig. 1. Experimental results for the amplitude of a bead at a single driving frequency: 
(a) shows the amplitude (in units of g) as a function of time, (b) is the power 
spectra for this time trace. The bead is a distance of 4 cm from the driving plate. 
The dotted line has slope of 2. (Data reproduced with permission of the authors 
of ref. |ig.) 

Fig. 2. Experimental results for the frequency response of a bead as a function of driving 
frequency. The bead is again a distance of 4 cm from the driving plate. The 
two curves result from successive measurements on an undisturbed system. The 
second curve is displaced downward so that the two are distinguishable. (Data 
reproduced with permission of the authors of ref. |10| .) 

Fig. 3. Basic ball and spring model for the system. The wall on the right oscillates to 
provide the driving force. 

Fig. 4. Frequency response of a the homogeneous system: (a) shows the response of the 
system with a damping coefficient (3 = 0, (b) shows the response with damping 
coefficient (3 = .003. 

Fig. 5. The 20 x 20 ball and spring system used for these simulations. The intensity 
indicates the strength of the bond connecting two balls: black indicates that 
Kh = K max , white indicates Kh = 0, with a linear gray scale for intermediate 
values. The large black ball at (8, 10) indicates the location where all of the 
measurements occur. 

Fig. 6. Frequency response of ball at (8, 10) as a function of frequency: (a) for low 
frequencies, (b) for intermediate frequencies and (c) for high frequencies. The 
pluses along the x-axis indicate the frequencies of the normal modes. 

Fig. 7. Phase of a single ball for all driving frequencies. There is a linear regime in the 
intermediate frequency range. The pluses along the x-axis mark the frequencies 
of the normal modes. 

Fig. 8. Typical modes for the system in three frequency ranges: (a) u = 0.1013c<; max ,; 

the mode is extended; (b) u = .6053c<; rnax ; the mode shows some localization; (c) 
u = 1.000w max ; the mode is very localized. 

Fig. 9. Values for the |cj|'s as a function of the mode frequency uji. The driving frequency 
is uj = .6uj max . 

Fig. 10. The most important modes. This figure shows the value of uii versus the driving 
frequency uj for the mode which has the largest value of |cj| at each value of value 



for uj. Note that above uj = uj max , the linear behavior starts to break down and 
the curve becomes very irregular. 

Fig. 11. The value of \ci\ as a function of the driving frequency uj for three successive 
modes of the system. 

Fig. 12. Results for the time trace of the amplitude as the bond network is changed: (a) 
is a sample of the time trace data, (b) shows the corresponding power spectra. A 
line with slope 2 is shown for comparison. 



